/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 Copyright (C) 2008, Proverbio Alessandro and Radice David

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

/** \file poly.cpp
    \brief Polynomial class and companion functions definitions
*/

#include<cassert>
#include<cmath>

#include "poly.hpp"

namespace spectral {

// =======================================================
// Constructors and destructors
// =======================================================

Poly::Poly(): _P_degree(-1),
              _P_a(NULL) {}

Poly::Poly(int const n): _P_degree(n) {
    _P_a = new double[_P_degree+1];
    for (int i=0;i<_P_degree;++i) {
        _P_a[i]=0;
    }
}

Poly::Poly(double const & x): _P_degree(0) {
    _P_a = new double[1];
    *_P_a=x;
}

Poly::Poly(Poly const &p): _P_degree(p._P_degree) {
    _P_a=new double[_P_degree+1];
    for(int i=0;i<=_P_degree;++i){
        _P_a[i]=p._P_a[i];
    }
}

Poly & Poly::operator=(Poly const & p) {
    delete[] _P_a;
    _P_degree=p._P_degree;
    _P_a=new double[_P_degree+1];
    for(int i=0;i<=_P_degree;++i){
        _P_a[i]=p._P_a[i];
    }
    return *this;
}

Poly::Poly(double const * const & coefficients, int degree): _P_degree(degree) {
    _P_a=new double[_P_degree+1];
    for(int i=0;i<=_P_degree;++i){
        _P_a[i]=coefficients[i];
    }
}

Poly::~Poly() { delete[] _P_a; }

// =======================================================
// Polynomial pubblic interface
// =======================================================

inline int Poly::degree() const { return _P_degree; }

inline double spectral::Poly::a(int n) const {
    assert(n<=_P_degree && n>=0);
    return _P_a[n];
}

inline double & Poly::operator[](int n) const {
    assert(n<=_P_degree && n>=0);
    return _P_a[n];
}

void Poly::setdegree(int n) {
    delete[] _P_a;
    _P_degree=n;
    _P_a=new double[_P_degree+1];
    for(int i=0;i<_P_degree+1;++i) {
        _P_a[i]=0.0;
    }
}

void Poly::setdegree() {
    while( ( std::abs( _P_a[_P_degree] )< POLYEPS ) && (_P_degree>=0) ) {
        --_P_degree;
    }

    double * tmp;
    tmp=new double[_P_degree+1];
    for(int i=0;i<=_P_degree;++i) {
        tmp[i]=_P_a[i];
    }
    delete[] _P_a;
    _P_a=tmp;
}

// =======================================================
// Special Operations
// =======================================================

Poly Poly::differentiate() const {
    double * tmp;
    tmp=new double[_P_degree];
    for(int i=0;i<_P_degree;++i){
        tmp[i]=double(i+1)*_P_a[i+1];
    }
    Poly P(tmp,_P_degree-1);
    delete[] tmp;
    return P;
}

Poly Poly::primitive() const {
    double * tmp;
    tmp=new double[_P_degree+2];
    tmp[0]=0.0;
    for(int i=0;i<=_P_degree;++i){
        tmp[i+1]=1.0/double(i+1)*_P_a[i];
    }
    Poly P(tmp,_P_degree+1);
    delete[] tmp;
    return P;
}

// This method will not call the primitive and evaluate methods
// for efficiency purpuses.
double Poly::integrate(double const & a, double const & b) const{
    double * tmp;
    tmp=new double[_P_degree+2];           //Evaluate the primitive
    tmp[0]=0.0;
    for(int i=0;i<=_P_degree;++i){
        tmp[i+1]=1.0/double(i+1)*_P_a[i];
    }
    double pb=tmp[_P_degree+1];
    double pa=tmp[_P_degree+1];
    double temp(0.0);
    for(int i=_P_degree;i>=0;--i){       //Horner algorithm
        temp=tmp[i];
        pb=pb*b+temp;
        pa=pa*a+temp;
    }
    delete[] tmp;
    return pb-pa;
}

// =======================================================
// Standard Polynomial Operations
// =======================================================

Poly & Poly::operator+=(Poly const &p){
    if(_P_degree<p._P_degree){
        double * coef;
        coef=new double[p._P_degree+1];
        for(int i=0;i<=_P_degree;++i){
            coef[i]=_P_a[i]+p._P_a[i];
        }
        for(int i=_P_degree+1;i<=p._P_degree;++i){
            coef[i]=p._P_a[i];
        }
        delete[] _P_a;
        _P_a=coef;
        _P_degree=p._P_degree;
    }
    else{
        for(int i=0;i<=p._P_degree;++i){
            _P_a[i]+=p._P_a[i];
        }
    }
    setdegree();
    return *this;
}

Poly & Poly::operator*=(double const &alpha){
    for(int i=0;i<=_P_degree;++i){
        _P_a[i]*=alpha;
    }
    return *this;
}

Poly & Poly::operator-=(Poly const & p){
    *this*=(-1);
    *this+=p;
    *this*=(-1);
    return *this;
}

bool Poly::operator==(Poly const & p) const {
    if(_P_degree!=p._P_degree) {
        return false;
    }
    for(int i=0;i<=_P_degree;++i) {
        if(std::abs(_P_a[i]-p._P_a[i]) > POLYEPS) {
            return false;
        }
    }
    return true;
}

// =======================================================
// External Functions
// =======================================================

Poly deflate(Poly const &p, double const & r) {
    int n=p.degree();
    Poly q(n-1);
    q[n-1]=p.a(n);
    for (int i=n-2;i>=0;--i) {
        q[i]=p.a(i+1)+q.a(i+1)*r;
    }
    return q;
}

Poly operator-(Poly p) {
    return Poly(double(0))-p;
}

Poly operator+(Poly a, Poly b) {
    return a+=b;
}

Poly operator-(Poly a, Poly b) {
    return a-=b;
}

Poly operator*(Poly a, double x) {
    return a*=x;
}

Poly operator*(double x, Poly a) {
    return a*=x;
}

double newton(Poly const & p,
              double delta,
              double eps,
              double x0,
              int nmax) {

    Poly dp=p.differentiate();
    double dlt=1.0;
    double xold;
    double xnew(x0);
    double res=std::abs(p(xnew));
    int k(1);

    while ( res>eps && dlt>delta && k<nmax ) {
        xold=xnew;
        assert(std::abs(dp(xold))>1e-16);
        xnew=xold-p(xold)/dp(xold);
        res=std::abs(p(xnew));
        dlt=std::abs(xnew-xold);
        ++k;
    }

    return xnew;
}

double * horner(Poly const & p,
            double delta,
            double eps,
            int nmax) {

    double * roots;
    int n=p.degree();
    roots=new double[n];
    double r;
    Poly q(p);

    for(int idx=0;idx<n;++idx) {
        r=newton(q,1e3*delta,1e3*eps,0.01,nmax);
        r=newton(p,delta,eps,r,nmax);
        q=deflate(q,r);
        roots[idx]=r;
    }
    return roots;
}

void legendre(Poly & L) {
    int n=L.degree();
    assert( n>0 );
    double * P=new double[n];
    double * G=new double[n];
    for (int i=0;i<n-1;++i) {
    P[i]=0;
    G[i]=0;
    }
    G[n-1]=0;
    for (int i=0;i<n+1;++i) {
        L[i]=0;
    }
    P[n-1]=0;
    L[1]=1;
    P[0]=1;
    for(int i=2;i<n+1;++i) {
    for(int j=0;j<i-1;++j) {
        G[j]=P[j];
    }
    for(int j=0;j<i;++j) {
        P[j]=L[j];
    }
    int k=0;
    while(i-k>0) {
        L[i-k]=(2*i-1)*P[i-k-1]/i;
        ++k;
    }
    int f=0;
    L[0]=-G[0]*(i-1)/i;
    while(i-f-1>0) {
        L[i-f-1]=L[i-f-1]-G[i-f-1]*(i-1)/i;
        ++f;
    }
    }
    delete[] P;
    delete[] G;
}

std::ostream & operator<<(std::ostream & strm,Poly P) {
    for(int i=P._P_degree; i>0; --i) {
        if( i==P._P_degree && P[i]<0 ) {
            strm << "- " << std::abs(P[i]) << " x^" << i;
        }
        else {
            strm << std::abs(P[i]) << " x^" << i;
        }
        if(P[i-1]>=0) {
            strm << " + ";
        }
        else {
            strm << " - ";
        }
    }
    strm << std::abs(P[0]);
    return strm;
}

std::istream & operator>>(std::istream & strm,Poly & P) {
    for(int i=0; i<P._P_degree+1; ++i) {
    strm >> P[i];
    }
    return strm;
}

} //namespace
